function [ results_fnames ] = osl_run_first_level_epoched( oat )

% [ results_fnames ] = osl_run_first_level_epoched( oat )
%
% takes in an OAT, which needs to be setup by calling oat=osl_setup_oat(S), struct
% and runs first level trialwise GLM
%
% This function should normally be called using osl_run_oat(oat);
%
% MWW 2012

global OSLDIR;

do_glm_statewise=oat.first_level.hmm_do_glm_statewise;

first_level=oat.first_level;
source_recon_name='source_recon';

if strcmp(oat.source_recon.modalities{1},'EEG')   % added by DM
    modality_meeg='EEG';
else
    modality_meeg='MEG';
end

% check if user has asked not to do GLM
if(~isfield(first_level,'doGLM')),
    first_level.doGLM = 1;
end

if(~isfield(first_level,'bc_trialwise')),
    first_level.bc_trialwise=0;
end;

if(~isfield(first_level, 'trialwise_directory')),
    first_level.trialwise_directory = [];
end

if isfield(first_level,'time_downsample_factor') && ~isempty(first_level,'time_downsample_factor'),
    error('first_level.time_downsample_factor no longer supported');
end;

if isfield(first_level,'tf_downsample_factor') && ~isempty(first_level,'tf_downsample_factor'),
    error('first_level.tf_downsample_factor no longer supported');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Set the downsample factors
% tf_downsample_factor is used in source space after recon weights are applied
% time_downsample_factor is used after TF transform (as part of osl_tf_transform call) in sensor space on complex values
if ~isempty(first_level.post_tf_downsample_factor)
    post_tf_ds_factor = 1./first_level.post_tf_downsample_factor;
else
    post_tf_ds_factor = 1;
end

if ~isempty(first_level.post_movingaverage_downsample_factor)
    post_movingaverage_ds_factor = 1./first_level.post_movingaverage_downsample_factor;
else
    post_movingaverage_ds_factor = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% set contrasts

contrast_list=first_level.contrast;
for c=1:length(contrast_list),
    contrast_list{c}=contrast_list{c}(:);
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check some of the settings

if(isfield(first_level,'connectivity_seed_mni_coord') && ~isempty(first_level.connectivity_seed_mni_coord)),
    error('This func does not handle seed connectivity analysis');
end;

if(sum(first_level.bc)>0 && first_level.baseline_timespan(1)>first_level.baseline_timespan(2)),
    error('first_level.baseline_timespan(1)>first_level.baseline_timespan(2)');
end;

if(~isfield(oat.first_level,'design_matrix') && ~isfield(oat.first_level,'design_matrix_summary')),
    error('Design matrix is not specified');
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% loop over sessions

clear stats;

for subi_todo=1:length(first_level.sessions_to_do),
    
    sub=first_level.sessions_to_do(subi_todo);
    
    %sub/length(first_level.submatfiles),
    
    %% load in the beamformer result for this session
    source_recon_results_fname = oat.source_recon.results_fnames{sub};
    source_recon_results=osl_load_oat_results(oat,source_recon_results_fname);
    
    % get D and update path in case OAT dir has been moved
    D = osl_get_oat_sensordata( source_recon_results );
    
    %% setup things that are common to all sessions
    if(subi_todo==1),
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up mask indices from mask/coords
        if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
            
            disp('Working in sensor space');
            
            chanindmeg = strmatch(modality_meeg, D.chantype);
            first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
            
        else,
            
            if(isfield(first_level,'mni_coords'))
                
                current_level_mni_coord = first_level.mni_coords;
                
                % load in lower level coords
                lower_level_mni_coord=source_recon_results.mni_coord;
                %lower_level_mask_fname  = [oat.source_recon.dirname '/' source_recon_name '_mask'];
                %[ lower_level_mni_coord, ~ ] = osl_mnimask2mnicoords(lower_level_mask_fname);
                
                [~, iA, mask_indices]=intersect(current_level_mni_coord, lower_level_mni_coord,'rows');
                
                % sort indices to be in the order of the lower level mni_coords:
                [~, gg]=sort(iA);
                indices_in_lower_level=mask_indices(gg);
                
                first_level_results.mask_indices_in_source_recon = indices_in_lower_level;
                first_level_results.mni_coord = lower_level_mni_coord(indices_in_lower_level,:);
                
                % the mni coordinates from the specified list
                str=['Desired MNI coordinates ' ];
                for vox=1:size(current_level_mni_coord,1),
                    str=[str ', [' num2str(current_level_mni_coord(vox,:)) ']'];
                end;
                disp(str);
                
                % the mni coordinates from the lower level mask
                str=['Nearest dipoles from lower level at MNI coordinates ' ];
                for vox=1:length(indices_in_lower_level),
                    str=[str ', [' num2str(lower_level_mni_coord(indices_in_lower_level(vox),:)) ']'];
                end;
                disp(str);
                
            else
                % setup std space brain
                first_level_results.gridstep=source_recon_results.gridstep;
                
                S=[];
                S.lower_level_mask_fname=[oat.source_recon.dirname '/' source_recon_name '_mask'];
                S.current_level_mask_fname=[oat.source_recon.dirname '/' oat.first_level.name '_mask'];
                S.current_level=first_level;
                S.lower_level_mni_coord=source_recon_results.mni_coord;
                S.lower_level_gridstep=source_recon_results.gridstep;
                
                [first_level_results.mask_indices_in_source_recon, first_level_results.mni_coord]=setup_mask_indices(S);
                lower_level_mni_coord=source_recon_results.mni_coord;
                
                clear S;
            end;
            
        end;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up potential time frequency bases in tf_settings, note that any time or
        % frequency averaging will be done AFTER the time-frequency decomposition
        
        if(~isempty(first_level.tf_hilbert_freq_ranges)),
            disp('oat.first_level.tf_hilbert_freq_ranges is set, freq_range, tf_freq_res and tf_num_freqs will be overridden');
            
            first_level.tf_num_freqs=size(first_level.tf_hilbert_freq_ranges,1);
            freq_range=[];
        else
            freq_range = first_level.tf_freq_range; % range used for stat testing
            if isempty(freq_range)
                % if not set as a first_level option, default to the range used
                % for source_recon
                freq_range=source_recon_results.source_recon.freq_range;
                
                if isempty(freq_range)
                    freq_range=[1 D.fsample/3];
                    disp(['WARNING: No frequency range has been specified. Set oat.source_recon.freq_range or oat.first_level.tf_freq_range. Using ' num2str(freq_range) ' Hz.']);
                end % isempty(freq_range)
            end % if isempty(freq_range)
            
        end;
        
        if ~isempty(freq_range),
            if(freq_range(2)>post_tf_ds_factor*D.fsample/2),
                warning(['Nyquist=' num2str(post_tf_ds_factor*D.fsample/2) ', freq_range(2)>post_tf_ds_factor*D.fsample/2. Results will be dodgey due to looking at frequencies higher than the Nyquist limit.']);
            end;
        else
            if(max(squash(first_level.tf_hilbert_freq_ranges))>post_tf_ds_factor*D.fsample/2),
                warning(['Nyquist=' num2str(post_tf_ds_factor*D.fsample/2) ', freq_range(2)>post_tf_ds_factor*D.fsample/2. Results will be dodgey due to looking at frequencies higher than the Nyquist limit.']);
            end;
        end;
        
        
        D_time_indices=find(source_recon_results.samples2use); % time indices in D.time used by source recon
        D_times = D.time(D_time_indices);
        
        if isempty(first_level.time_range)
            first_level.time_range=source_recon_results.woi;
        end
        
        tf_settings = [];
        tf_settings.tf_method               = first_level.tf_method;
        tf_settings.tf_logtransform         = first_level.tf_logtransform;
        tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
        tf_settings.tf_hilbert_freq_ranges  = first_level.tf_hilbert_freq_ranges;
        tf_settings.tf_freq_range           = freq_range; % set above
        tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
        tf_settings.tf_multitaper_ncycles   = first_level.tf_multitaper_ncycles;
        tf_settings.raw_times               = D_times;
        tf_settings.ds_factor               = post_tf_ds_factor;
        tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
        tf_settings.tf_multitaper_taper     = first_level.tf_multitaper_taper;
        tf_settings.tf_multitaper_twin      = first_level.tf_multitaper_twin;
        tf_settings.tf_multitaper_freqsmooth= first_level.tf_multitaper_freqsmooth;
        tf_settings.tf_calc_amplitude=0; % we will calc the amplitude after applying recon weights
        tf_settings.tf_hilbert_do_bandpass_for_single_freq = first_level.tf_hilbert_do_bandpass_for_single_freq;
        
        out = osl_tf_transform( tf_settings ); % passed without a data matrix, osl_tf_transform just gets container sizes
        
        if isfield(out, 'tf_morlet_basis')
            first_level.tf_morlet_basis = out.tf_morlet_basis;
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup time and frequency windows/indices for using later
        
        % setup result time dimension
        %tf_time_indices_into_D_times = ( out.tf_times >= first_level.time_range(1) ) & ( out.tf_times <= first_level.time_range(2) );
        tf_time_indices_into_D_times = false(1,numel(out.tf_times));
        for i=1:size(first_level.time_range,1)
            tf_time_indices_into_D_times(( out.tf_times >= first_level.time_range(i,1) ) & ( out.tf_times <= first_level.time_range(i,2) )) = true;
        end
        
        first_level_results.times = out.tf_times(tf_time_indices_into_D_times);
        % downsample first_level_results.times
        % temporally downsample
        if post_movingaverage_ds_factor~=1,
            tmp=ones(size(first_level_results.times));
            tmp=osl_spm_resample(tmp,post_movingaverage_ds_factor);
            first_level_results.times=linspace(first_level_results.times(1),first_level_results.times(end),length(tmp));
        end;
        
        if first_level.time_average
            first_level_results.times=mean(first_level_results.times);
            tres=1;
        end
        
        % setup result freq dimension
        first_level_results.frequencies = out.tf_freqs;
        if first_level.freq_average,
            first_level_results.frequencies=mean(first_level_results.frequencies);
        end
        
        % setup result spatial dimension
        Nvoxels_out=length(first_level_results.mask_indices_in_source_recon);
        
        clear out
        
        %%%%%%%%%%%%%%%%
        %% results containers
        ntpts                       = length(first_level_results.times);
        nfreqs                      = length(first_level_results.frequencies);
        
        artefact_map = nan(Nvoxels_out,nfreqs);
        
        if(first_level.doGLM)
            first_level_results.stdcope=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs]);% num_voxels x  num_timepoints x num_contrasts x num_freqs
            first_level_results.cope=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs]);% num_voxels x num_timepoints x num_contrasts x num_freqs
            
        end;
        
        first_level_results.pseudo_zstat_var = zeros([Nvoxels_out,1],'single');
        first_level_results.wnorm_nai = cell(1,Nvoxels_out);
        
    end; % end of if(sub==1)
    
    disp(['Reconstruct time courses and computing stats for dataset ' source_recon_results.fname]);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% setup HMM
    if strcmp(source_recon_results.recon_method,'none'),
        NKrecon=1;
        NKglm=NKrecon;
    else
        NKrecon=numel(source_recon_results.BF.inverse.W.MEG);
        NKglm=NKrecon;
        
        if(~do_glm_statewise),
            NKglm=1;
        end;
        
        if(NKglm>1)
            if post_movingaverage_ds_factor~=1,
                error('Can not do any downsampling when using multiple classes (e.g. HMM)');
            end;
        end;
    end;
    
    if(first_level.doGLM && NKglm>1)
        first_level_results.stdcope_by_state=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs,NKglm]);% num_voxels x  num_timepoints x num_contrasts x num_freqs x num_states
        first_level_results.cope_by_state=zeros([Nvoxels_out,ntpts,length(first_level.contrast),nfreqs,NKglm]);% num_voxels x num_timepoints x num_contrasts x num_freqs x num_states
    end;
    
    %%%%%%%%%%%%%%%
    %% setup channels
    chanindmeg = strmatch(modality_meeg, D.chantype);
    chanind = setdiff(chanindmeg, D.badchannels);
    if isempty(chanind)
        error(['No good ' modality_meeg ' channels were found.']);
    end
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
        % for sensor space - also include bad sensors sensors
        first_level_results.mask_indices_in_source_recon=1:length(chanindmeg);
        chanind = chanindmeg;
    end;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% setup the GLM design matrix for this session
    % setup trial types
    trialtypes=[];
    triallist=[];
    
    for i=1:length(source_recon_results.source_recon.conditions), % indexes conditions/triggers within subset
        
        trigname=source_recon_results.source_recon.conditions{i};
        Ntrialspercond=length(D.pickconditions(trigname)); %% number of trials for this condition
        triallist=[triallist , D.pickconditions(trigname)];
        trialtypes=[trialtypes ; ones(Ntrialspercond,1).*i];
        
    end;
    
    S3=[];
    S3.trialtypes=trialtypes;
    
    if numel(first_level.design_matrix_summary)==length(oat.source_recon.results_fnames) &&(ischar(first_level.design_matrix_summary{sub})),
        
        S3.Xsummary=first_level.design_matrix_summary{sub};
        
        if(isfield(first_level,'trial_rejects'))
            S3.trial_rejects=first_level.trial_rejects;
        else
            warning('oat.first_level.trial_rejects not set');
        end;
    else
        S3.Xsummary=first_level.design_matrix_summary;
    end;
    
    first_level.x=setup_beamformer_designmatrix(S3);
    
    const_reg=0;
    for p=1:size(first_level.x,2),
        if(range(first_level.x(:,p))==0),
            if(const_reg),
                warning('Multiple constant regressors in the design matrix!!!!?');
            end;
            const_reg=1;
            
        end;
        
        if(~any(first_level.x(:,p))),
            warning(['Empty (all zeros) regressor (EV' num2str(p) ') in the design matrix!!!!? Could cause problems for subsequent OLS group analysis.']);
        end;
    end;
    
    if(first_level.do_glm_demean)
        warning('CAREFUL: are you sure you want the first_level.do_glm_demean flag off? This means that main effect contrasts will not be interpetable as mean effects');
    end;
    
    if(const_reg && first_level.do_glm_demean)
        warning('CAREFUL: are you sure you want the first_level.do_glm_demean flag on with constant regressors in the design matrix!!!!?');
    elseif(~const_reg && ~first_level.do_glm_demean)
        warning('CAREFUL: are you sure you want no first_level.do_glm_demean flag on with no constant regressors in the design matrix!!!!?');
    end;
    
    Ntrials   = length(triallist);
    
    %%%%%%%%%%%%%%%
    %% load in sensor data
    first_level_results.D_sensor_data=D;
    first_level_results.chanind=chanind;
    
    sensor_data=D(chanind, D_time_indices, triallist);
    
    classchanind=find(strcmp(D.chanlabels,'Class'));
    if(isempty(classchanind)),
        error(['No ''Class'' chanlabel in: ' D.fname]);
    end;
    
    % class_samples_inds_recon is used for using weights in recon
    for si=1:NKrecon,
        class_samples_inds_recon{si} = (D(classchanind, D_time_indices, triallist)==si);
    end;
    
    % class_samples_inds_glm is used for doing GLM state(class)-wise
    if(~do_glm_statewise),
        class_samples_inds_glm{1}=ones(1,length(D_time_indices),length(triallist));
    else
        class_samples_inds_glm=class_samples_inds_recon;
    end;
    
    for si=1:NKrecon,
        state_tinds=class_samples_inds_recon{si};
        
        total_time=(D_times(end)-D_times(1))*size(state_tinds,3);
        state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
        txt=['State ' num2str(si) ' is active for ' num2str(state_time) 'secs'];
        
        if(state_time<10),
            warning([txt '. It will be excluded from the GLM analysis.']);
        else
            disp(txt);
        end;
    end;
    clear state_tinds;
    
    %%%%%%%%%%%%%%%
    %% Do T-F transform in sensor space
    %% Transform the data - either downsampling in time domain, or a TF transform
    
    if(~strcmp(first_level.tf_method,'none'))
        disp('Doing T-F transform in sensor space...');
    end;
    
    tf_settings = [];
    tf_settings.tf_method               = first_level.tf_method;
    tf_settings.tf_logtransform         = first_level.tf_logtransform;
    tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
    tf_settings.tf_hilbert_freq_ranges  = first_level.tf_hilbert_freq_ranges;
    tf_settings.tf_freq_range           = freq_range; % set above
    tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
    tf_settings.raw_times               = D_times;
    tf_settings.ds_factor               = post_tf_ds_factor;
    tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
    tf_settings.tf_multitaper_taper     = first_level.tf_multitaper_taper;
    tf_settings.tf_multitaper_twin      = first_level.tf_multitaper_twin;
    tf_settings.tf_multitaper_freqsmooth = first_level.tf_multitaper_freqsmooth;
    tf_settings.tf_calc_amplitude=0; % we will calc the amplitude after applying recon weights
    tf_settings.tf_hilbert_do_bandpass_for_single_freq = first_level.tf_hilbert_do_bandpass_for_single_freq;
    
    if isfield(first_level,'tf_morlet_basis')
        tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
    end
    
    dims=size(sensor_data);
    dims_tf=[dims(1),length(tf_time_indices_into_D_times),Ntrials,first_level.tf_num_freqs];
    disp(['size(sens_data_tf)=' mat2str((dims_tf)) ', (nsens x ntpts x ntri x nfreq)']);
    
    if(prod(size(dims_tf)) > 1e9)
        warning(['size(sens_data_tf)=' mat2str((dims_tf)) '. May cause RAM problems. Consider using less freq bins or using oat.first_level.post_tf_downsample_factor>1.']);
    end;
    
    sensor_data_tf=zeros(dims_tf);
    
    if(first_level.tf_num_freqs>1), ft_progress('init', 'etf');end;
    for indind=1:size(sensor_data,1),
        if(first_level.tf_num_freqs>1), ft_progress(indind/size(sensor_data,1));end;
        
        for tri=1:size(sensor_data,3),
            out = osl_tf_transform( tf_settings , sensor_data(indind,:,tri) );
            sensor_data_tf(indind,:,tri,:)=permute(out.dattf,[2 3 1]);
        end;
    end;
    
    % get indices into the class vectors, in order to downsample as per
    % osl_tf_transform
    ds_inds = round(linspace(1,numel(D_times),numel(out.tf_times)));
    for si=1:NKglm,
        class_samples_inds_glm_tf{si}   = class_samples_inds_glm{si}(:,ds_inds,:);
        class_samples_inds_recon_tf{si} = class_samples_inds_recon{si}(:,ds_inds,:);
    end;
    
    if(first_level.tf_num_freqs>1), ft_progress('close');end;
    
    %%%%%%%%%%%%%%%%
    %% Now do windowing over time
    sensor_data_tf=sensor_data_tf(:,tf_time_indices_into_D_times,:,:);
    for si=1:NKglm,
        class_samples_inds_glm_tf{si}=class_samples_inds_glm_tf{si}(:,tf_time_indices_into_D_times,:);
    end;
    for si=1:NKrecon,
        class_samples_inds_recon_tf{si}=class_samples_inds_recon_tf{si}(:,tf_time_indices_into_D_times,:);
    end;
    
    %%%%%%%%%%%%%%%%
    %% Temporally downsample state time courses (data will be done later after taking abs (for t-f data) inside the voxel loop)
    for si=1:NKglm,
        
        state_tinds=class_samples_inds_glm_tf{si};
        
        % moving average within the window specified
        if ~isempty(first_level.time_moving_av_win_size)
            
            if(NKglm>1)
                error('first_level.time_moving_av_win_size not compatible with more than one class');
            end;
            
            % do nothing
        end
        
        % temporally downsample
        if post_movingaverage_ds_factor~=1,
            clear state_tinds_new;
            for tri=1:size(state_tinds,3),
                tmp=permute(state_tinds(1,:,tri),[1,2,3]);
                tmp=osl_spm_resample(tmp,post_movingaverage_ds_factor);
                state_tinds_new(1,:,tri)=tmp;
            end;
            state_tinds=state_tinds_new;
            
            % need binary states
            state_tinds=round(state_tinds);
            state_tinds(state_tinds<0)=0;
            state_tinds(state_tinds>1)=1;
        end;
        
        class_samples_inds_glm_tf{si}=state_tinds;
    end;
    
    % no need to temporally downsample or moving average class_samples_inds_recon_tf
    
    % downsample first_level_results.times
    
    
    %%%%%%%%%%%%%%%%
    %% setup design matrices that can be reused over voxels
    ntpts=length(first_level_results.times);
    x_stored=cell(NKglm,ntpts);
    pinvxtx_stored=cell(NKglm,ntpts);
    pinvx_stored=cell(NKglm,ntpts);
    contrast_list_mat=cell2mat(contrast_list);
    contrast_list_mat2=eye(size(first_level.x,2));
    
    isrepresented = zeros(NKglm,size(state_tinds,2));
    for si=1:NKglm,
        
        state_tinds=class_samples_inds_glm_tf{si};
        
        total_time=(D_times(end)-D_times(1))*size(state_tinds,3);
        state_time=total_time*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
        
        if(state_time>10),
            xc=[];
            contrast_list_matc=[];
            contrast_list_matc2=[];
            
            
            for t = 1 : size(state_tinds,2), % indexes time points within trials
                
                %t=nearest(D_times,first_level_results.times(ts));
                
                %disp(num2str(t/ntpts));
                
                % x is (ntrials x nevs)
                x=first_level.x;
                
                % check if the state is associated with any data for this
                % timepoint
                if isempty(find(permute(state_tinds(1,t,:),[3 1 2])))
                else
                    isrepresented(si,t) = 1;
                end
                
                x=x(find(permute(state_tinds(1,t,:),[3 1 2])),:);
                
                if(oat.first_level.do_glm_demean)
                    for pp=1:size(x,2),
                        x(:,pp)=x(:,pp)-mean(x(:,pp));
                    end;
                end;
                
                x_stored{si,t}=x;
                
                pinvxtx_stored{si,t}=pinv(x'*x);
                pinvx_stored{si,t}=pinv(x);
                
                xc=sparse(blkdiag(xc,x));
                
                contrast_list_matc=sparse(blkdiag(contrast_list_matc,contrast_list_mat));
                
                contrast_list_matc2=sparse(blkdiag(contrast_list_matc2,contrast_list_mat2));
                
            end;
            
            
            contrast_list_matc=sparse(contrast_list_matc);
            contrast_list_matc2=sparse(contrast_list_matc2);
            
            % xc is (blkdiag of x from repeated for each time point)
            % x is ([ntpts*ntrials] x [ntpts*nevs])
            
            xc_stored{si}=sparse(xc);
            xcxc=full(xc'*xc);
            pinvxctxc_stored{si}=sparse(pinv(xcxc));
            xc=full(xc);
            pinvxc_stored{si}=sparse(pinv(xc));
        end;
    end;
    clear state_tinds;
    
    % check that the non-artefact states account for all timepoints
    if ~isempty(oat.source_recon.artefact_chanlabel)
        iscaptured = sum(isrepresented(1:size(isrepresented,1)-1,:),1);
        if any(~iscaptured)
            percentmissing = 100 * (sum(~iscaptured) / numel(iscaptured));
            warning(['The non-artefact state(s) fail to cover ' num2str(percentmissing) '% of timepoints within the trials']);
        end
    end
    
    %%%%%%%%%%%%%%%%
    %% do everything from now on one voxel at a time (for the sake of
    %% RAM) in the first level mask space
    if(first_level.doGLM)
        disp('First level COPEs outputted will have dimension Nvoxels x Ntpts x Ncontrasts x Nfreqs:');
        tmp=[size(first_level_results.cope) 1 1];
        disp([num2str(tmp(1:4))]);
    end;
    
    % containers
    if(first_level.doGLM)
        covarpe_by_state=zeros([ntpts,size(first_level.x,2),size(first_level.x,2),length(first_level_results.frequencies),NKglm]);% num_timepoints x num_contrasts x num_contrasts x num_freqs x num_states
        pe_by_state=zeros([ntpts,size(first_level.x,2),length(first_level_results.frequencies),NKglm]);% num_timepoints x num_contrasts x num_freqs x num_states
    end;
    
    if first_level.save_trialwise_data
        % transform-specific parameters
        switch first_level.tf_method
            case 'none'
                trlwise.type = 'time_domain';
                trlwise.dat     = nan(Nvoxels_out,Ntrials,ntpts,nfreqs);
            otherwise
                trlwise.type    = first_level.tf_method;
                trlwise.dat     = nan(Nvoxels_out,Ntrials,ntpts,nfreqs);
                trlwise.freqs   = first_level_results.frequencies;
                trlwise.times   = first_level_results.times;
        end
        trlwise.triallist = triallist;
        trlwise.rawdat    = fullfile(D.path,D.fname);
        % spatial parameters
        if (isfield(first_level,'mni_coords'))
            trlwise.mni_coords   = first_level_results.mni_coord;
        else
            trlwise.mask         = read_avw([oat.source_recon.dirname '/' oat.first_level.name '_mask']);
            trlwise.gridstep     = first_level_results.gridstep;
        end
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% now loop over voxels
    ft_progress('init', 'etf');
    
    % debug:
    %[indfind, mni_coord]=nearest_vec(source_recon_results.mni_coord,[-32 -26 50]);
    
    S2=[];
    S2.source_recon_results=source_recon_results;
    S2.chanind=chanind;
    S2.class_samples_inds=class_samples_inds_recon_tf;
    S2.num_freqs=first_level.tf_num_freqs;
    %S2.sensor_data_tf=sensor_data_tf;
    
    Ntrials=size(class_samples_inds_recon{1},3);
    for kk=1:NKrecon,
        for tri=1:Ntrials,
            S2.timeinds{kk,tri}=find(class_samples_inds_recon_tf{kk}(1,:, tri)); % time indices for class kk
            S2.sensor_data_sub{kk,tri}=reshape(sensor_data_tf(:,S2.timeinds{kk,tri},tri,:),[size(sensor_data_tf,1),length(S2.timeinds{kk,tri})*size(sensor_data_tf,4)]);
        end;
    end;
    
    %clear sensor_data;
    
    for indind=1:Nvoxels_out, % indexes brain space
        
        ft_progress(indind/Nvoxels_out);
        
        S2.voxind=first_level_results.mask_indices_in_source_recon(indind);
        
        [ dat_tf wnorms wnorm wnorms_nai wnorm_nai ] = get_voxel_recon_timecourse_vector( S2 );
        
        %dat_tf=dat_tf./sqrt(wnorms_nai);
        
        first_level_results.pseudo_zstat_var(indind)=var(squash(dat_tf));
        
        first_level_results.wnorm_nai{indind}=wnorm_nai;
        
        if sum(isnan(squash(dat_tf)))~=0,
            
            if(first_level.doGLM)
                first_level_results.cope(indind,:,:,:)=0;
                first_level_results.stdcope(indind,:,:,:)=inf;
                if(NKglm>1)
                    first_level_results.cope_by_state(indind,:,:,:,:)=0;
                    first_level_results.stdcope_by_state(indind,:,:,:,:)=inf;
                end;
            end;
            error('NANs in recon data');
            
        else
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %% calc the amplitude (up to now dat_tf is complex)
            if ~strcmp(first_level.tf_method,'none'),
                dat_tf=abs(dat_tf);
            end;
            
            % perform weights normalization
            dat_tf=dat_tf./sqrt(wnorms_nai);
            
            dat=dat_tf;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %% moving average within
            %% the window specified
            if ~isempty(first_level.time_moving_av_win_size)
                
                fsample=1/(D_times(2)-D_times(1));
                windowSize = first_level.time_moving_av_win_size*fsample;
                
                clear dat_new;
                for tri=1:size(dat,1),
                    for f=1:size(dat,3),
                        tmp=permute(dat(tri,:,f),[1,2,3]);
                        tmp=moving(tmp,round(windowSize));
                        dat_new(tri,:,f)=tmp;
                    end;
                end;
                dat=dat_new;
            end
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %% temporally downsample (needs to be done after abs (for t-f data))
            if post_movingaverage_ds_factor~=1,
                clear dat_new;
                for tri=1:size(dat,1),
                    for f=1:size(dat,3),
                        tmp=permute(dat(tri,:,f),[1,2,3]);
                        tmp=osl_spm_resample(tmp,post_movingaverage_ds_factor);
                        dat_new(tri,:,f)=tmp;
                    end;
                end;
                dat=dat_new;
                
            end;
            
            glm_input_times=linspace(first_level.time_range(1),first_level.time_range(2),length(dat));
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %% now do averaging over time and frequency
            if(first_level.time_average)
                dat=mean(dat,2);
            end;
            
            if(first_level.freq_average)
                dat=mean(dat,3);
            end;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if first_level.save_trialwise_data
                trlwise.dat(indind,:,:,:) = dat;
            end % if first_level.save_data
            
            if(first_level.bc_trialwise),
                baseline_time_indices=first_level_results.times<first_level.baseline_timespan(2) & first_level_results.times>first_level.baseline_timespan(1);
                
                if(sum(baseline_time_indices)>0)
                    % normalise by noise in baseline window:
                    for tri=1:size(dat,1), % indexes trials
                        for f = 1 : size(dat,3),
                            dat(tri,:,f)=(dat(tri,:,f)-mean(dat(tri,baseline_time_indices,f)));
                        end;
                    end;
                else
                    warning('Could not baseline correct as there are no timepoints in the specified baseline period');
                end;
            end;
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %% Do GLM
            if first_level.doGLM,
                npes=size(contrast_list_mat2,2);
                
                if(size(dat,1)>1)
                    for f = 1 : size(dat,3),
                        clear total_time state_time;
                        for si=1:NKglm,
                            
                            state_tinds=class_samples_inds_glm_tf{si};
                            
                            total_time{si}=(D_times(end)-D_times(1))*size(state_tinds,3);
                            state_time{si}=total_time{si}*sum(sum(squeeze(state_tinds)))/prod(size(squeeze(state_tinds)));
                            
                            if(state_time{si}<10),
                                
                                if(NKglm>1)
                                    first_level_results.cope_by_state(indind,:,:,f,si)=0;
                                    first_level_results.stdcope_by_state(indind,:,:,f,si)=inf;
                                end;
                            else
                                
                                dattfc=[];
                                for t = 1 : size(dat,2), % indexes time points within trials
                                    
                                    dattf=dat(find(permute(state_tinds(1,t,:),[3 1 2])),t,f);
                                    
                                    if(oat.first_level.do_glm_demean)
                                        dattf=dattf-mean(dattf);
                                    end;
                                    
                                    dattfc=cat(1,dattfc,dattf);
                                    
                                end;
                                
                                % xc is (blkdiag of x from repeated for each time point)
                                % x is ([ntpts*ntrials] x [ntpts*nevs])
                                x=xc_stored{si};
                                pinvxtx=pinvxctxc_stored{si};
                                pinvx=pinvxc_stored{si};
                                
                                [peout, covarpeout]=glm_fast_for_meg(dattfc,x,pinvxtx,pinvx);
                                
                                pe_by_state(:,:,f,si)=reshape(peout,npes,size(dat,2))';
                                
                                for t= 1 : size(dat,2),
                                    covarpe_by_state(t,:,:,f,si)=covarpeout((t-1)*npes+1:t*npes,(t-1)*npes+1:t*npes);
                                end;
                                
                            end;
                        end; % si (state index)
                        
                        %% do fixed effects averaging of states
                        % time loop
                        
                        
                        for iTime = 1 : size(pe_by_state,1), % indexes time points within trials
                            
                            if(NKglm>1)
                                pe_tmp=[];
                                covpe_tmp=[];
                                z=[];
                                
                                % one state may be an artefact state
                                if ~isempty(oat.source_recon.artefact_chanlabel)
                                    % the state with highest index is the
                                    % artefact state
                                    nstatess = size(pe_by_state,4) - 1; % don't use the artefact state
                                else
                                    nstatess = size(pe_by_state,4);
                                end
                                
                                for si = 1 : nstatess, % indexes state
                                    if(state_time{si}>=10),
                                        z=[z; eye(npes)];
                                        tmppe=permute(pe_by_state(iTime,:,f,si),[2 3 1]);
                                        pe_tmp=[pe_tmp; tmppe];
                                        
                                        if(strcmp(oat.first_level.cope_type,'coape')),
                                            tmppe=abs(tmppe);
                                        end;
                                        
                                        tmpcovpe=permute(covarpe_by_state(iTime,:,:,f,si),[2 3 4 1]);
                                        covpe_tmp=blkdiag(covpe_tmp,tmpcovpe);
                                        if(NKglm>1)
                                            first_level_results.cope_by_state(indind,iTime,:,f,si)=contrast_list_mat'*tmppe;
                                            first_level_results.stdcope_by_state(indind,iTime,:,f,si)=sqrt(diag(contrast_list_mat'*tmpcovpe*contrast_list_mat));
                                        end;
                                    end;
                                end;
                                
                                [gam, covgam] = fixed_effects(pe_tmp,z,covpe_tmp);
                                %gam=mean(cope_tmp);
                            else
                                gam=permute(pe_by_state(iTime,:,f),[2 1]);
                                covgam=permute(covarpe_by_state(iTime,:,:,f),[2 3 1]);
                            end;
                            
                            
                            if(NKglm>1)
                                % if there is enough data,
                                % average the artefact state over all
                                % timepoints and conditions
                                if ~isempty(oat.source_recon.artefact_chanlabel)
                                    if(state_time{nstatess+1}>=10)
                                        tmp = mean(pe_by_state(:,:,:,nstatess+1),2);
                                        artefact_map(indind,:) = squeeze(mean(tmp,1));
                                    else
                                    end
                                end
                            end
                            
                            if(strcmp(oat.first_level.cope_type,'coape')),
                                gam=abs(gam);
                            end;
                            
                            % NOW: apply contrasts
                            first_level_results.cope(indind,iTime,:,f)=contrast_list_mat'*gam;
                            
                            if(strcmp(oat.first_level.cope_type,'acope')),
                                first_level_results.cope(indind,iTime,:,f)=abs(first_level_results.cope(indind,iTime,:,f));
                            end;
                            
                            first_level_results.stdcope(indind,iTime,:,f)=sqrt(diag(contrast_list_mat'*covgam*contrast_list_mat));
                            
                            if(NKglm==1)
                                first_level_results.cope_by_state(indind,iTime,:,f,1)=first_level_results.cope(indind,iTime,:,f);
                                first_level_results.stdcope_by_state(indind,iTime,:,f,1)=first_level_results.stdcope(indind,iTime,:,f);
                            end;
                            
                        end; % for iTime
                        
                        %% do bc on the copes
                        if(sum(first_level.bc)>0),
                            baseline_time_indices=first_level_results.times<first_level.baseline_timespan(2) & first_level_results.times>first_level.baseline_timespan(1);
                            
                            for c=1:length(contrast_list),
                                if(first_level.bc(c))
                                    first_level_results.cope(indind,:,c,f)=first_level_results.cope(indind,:,c,f)-mean(first_level_results.cope(indind,baseline_time_indices,c,f),2);
                                end;
                            end;
                            
                        end;
                    end; % end freq for loop
                    
                    if(0)
                        indind
                        iContrast=[2,4,5,6];iFreq=2;
                        figure(1);iFreq=2;plot(squeeze(first_level_results.cope(indind,:,iContrast,iFreq)))
                        
                    end;
                else
                    first_level_results.cope(indind,:,1,:)=dat(:,:);
                    first_level_results.stdcope(indind,:,1,:)=1;
                    if(NKglm>1)
                        first_level_results.cope_by_state(indind,:,1,:)=dat(:,:);
                        first_level_results.stdcope_by_state(indind,:,1,:)=1;
                    end;
                end;
            else
                warning('GLM not run: first_level.doGLM == 0');
            end;
            
        end % if first_level.doGLM
    end % for indind = 1:Nvoxels_out
    
    if (first_level.doGLM && ~isempty(first_level.connectivity_seed_mni_coord) && NKglm>1),
        first_level_results.cope_by_state(first_level_results.seed_index,:,:,:,:)=0;
        first_level_results.stdcope_by_state(first_level_results.seed_index,:,:,:,:)=inf;
    end;
    
    ft_progress('close');
    
    
    %% combine the planar gradiometers if in sensorspace
    if strcmp(source_recon_results.recon_method,'none')
        switch first_level.sensor_space_combine_planars
            case 'combine_cartesian'
                
                % get labels for planar 1, planar 2, and magnetometers
                lay_fn = fullfile(OSLDIR,'layouts','neuromag306planar3.lay');
                cfg=[];
                cfg.layout = lay_fn;
                layplan1 = ft_prepare_layout(cfg);
                
                lay_fn = fullfile(OSLDIR,'layouts','neuromag306planar2.lay');
                cfg=[];
                cfg.layout = lay_fn;
                layplan2 = ft_prepare_layout(cfg);
                
                lay_fn = fullfile(OSLDIR,'layouts','neuromag306mag.lay');
                cfg=[];
                cfg.layout = lay_fn;
                laygrad = ft_prepare_layout(cfg);
                
                lay_fn = fullfile(OSLDIR,'layouts','neuromag306cmb.lay');
                cfg=[];
                cfg.layout = lay_fn;
                laycmb = ft_prepare_layout(cfg);
                
                % get indices for planar 1, planar 2, and magnetometers
                [planar1labels,planar1inds_inD] = intersect(D.chanlabels,layplan1.label);
                [planar2labels,planar2inds_inD] = intersect(D.chanlabels,layplan2.label);
                [maglabels,maginds_inD]         = intersect(D.chanlabels,laygrad.label);
                
                magcope             = first_level_results.cope(maginds_inD,:,:,:);
                magstd              = first_level_results.stdcope(maginds_inD,:,:,:);
                if(NKglm>1)
                    magstatecope        = first_level_results.cope_by_state(maginds_inD,:,:,:);
                    magstatestd         = first_level_results.stdcope_by_state(maginds_inD,:,:,:);
                end;
                magpseudozstatvar   = first_level_results.pseudo_zstat_var(maginds_inD,:,:,:);
                
                % adopting the order of the planar1's, step through the
                % planar2's and find the corresponding pair.  If not
                % present, throw the planar1 data away.  If present,
                % combine the data from the two channels
                pairlist                = [];
                combinedCope            = nan(numel(planar1labels),size(first_level_results.cope,2),size(first_level_results.cope,3),size(first_level_results.cope,4));
                %combinedCopeByState     = nan(numel(planar1labels),size(first_level_results.cope,2),size(first_level_results.cope,3),size(first_level_results.cope,4));
                combinedStdCope         = nan(numel(planar1labels),size(first_level_results.cope,2),size(first_level_results.cope,3),size(first_level_results.cope,4));
                %combinedStdCopeByState  = nan(numel(planar1labels),size(first_level_results.cope,2),size(first_level_results.cope,3),size(first_level_results.cope,4));
                count = 0;
                clear newcmblabels
                for iFindPair = 1:numel(planar1inds_inD)
                    correspondingInd = find(strncmpi(planar1labels{iFindPair},planar2labels,6));
                    if ~isempty(correspondingInd)
                        count = count + 1;
                        pairlist(iFindPair,:) = [iFindPair,correspondingInd];
                        
                        % combine the cope
                        planar1res_temp = first_level_results.cope(planar1inds_inD(iFindPair),:,:,:);
                        planar2res_temp = first_level_results.cope(planar2inds_inD(correspondingInd),:,:,:);
                        if (strcmp(oat.first_level.cope_type,'acope')) || (~(strcmp(oat.first_level.tf_method,'none'))) % if the copes are rectified
                            combinedCope(count,:,:,:) = (planar1res_temp+planar2res_temp)./2;
                            if count == 1
                                disp('Averaging the (A)COPEs for orthogonal planar gradiometers')
                            end
                        else
                            combinedCope(count,:,:,:) = sqrt( planar1res_temp.^2 + planar2res_temp.^2  );
                            if count == 1
                                warning('Gradiometers combined with ERF and COPE: gradiometer COPEs will be rectified');
                            end
                        end
                        
                        % combine the stdcope
                        planar1res_temp = first_level_results.stdcope(planar1inds_inD(iFindPair),:,:,:);
                        planar2res_temp = first_level_results.stdcope(planar2inds_inD(correspondingInd),:,:,:);
                        combinedStdCope(count,:,:,:) = sqrt( planar1res_temp.^2 + planar2res_temp.^2  ); % same transformation here but for different reasons
                        
                        % combine the cope_by_state
                        if(0),%if(NKglm>1)
                            planar1res_temp = first_level_results.cope_by_state(planar1inds_inD(iFindPair),:,:,:);
                            planar2res_temp = first_level_results.cope_by_state(planar2inds_inD(correspondingInd),:,:,:);
                            
                            if (strcmp(oat.first_level.cope_type,'acope')) || (~(strcmp(oat.first_level.tf_method,'none'))) % if the copes are rectified
                                combinedCopeByState(count,:,:,:) = (planar1res_temp+planar2res_temp)./2;
                                if count == 1
                                    disp('Averaging the (A)COPEs (by state) for orthogonal planar gradiometers')
                                end
                            else
                                combinedCopeByState(count,:,:,:) = sqrt( planar1res_temp.^2 + planar2res_temp.^2  );
                                %warning('Gradiometers combined with ERF and COPE: gradiometer COPEs will be rectified');
                            end
                            
                            % combine the stdcope_by_state
                            planar1res_temp = first_level_results.stdcope_by_state(planar1inds_inD(iFindPair),:,:,:);
                            planar2res_temp = first_level_results.stdcope_by_state(planar2inds_inD(correspondingInd),:,:,:);
                        end;
                        %combinedStdCopeByState(count,:,:,:) = sqrt( planar1res_temp.^2 + planar2res_temp.^2  );
                        
                        labind = find(strncmpi(planar1labels{iFindPair},laycmb.label,6));
                        newcmblabels{count} = laycmb.label{labind};
                    else
                        
                    end
                end % for iFindPair = 1:numel(planar1inds_inD)
                
                combinedCope(isnan(combinedCope)) = [];
                combinedStdCope(isnan(combinedStdCope)) = [];
                %combinedCopeByState(isnan(combinedCopeByState)) = [];
                %combinedStdCopeByState(isnan(combinedStdCopeByState)) = [];
                
                % re-baseline
                if(sum(first_level.bc)>0),
                    baseline_time_indices=first_level_results.times<first_level.baseline_timespan(2) & first_level_results.times>first_level.baseline_timespan(1);
                    
                    for c=1:length(contrast_list),
                        if(first_level.bc(c))
                            combinedCope(:,:,c,:) = combinedCope(:,:,c,:) - repmat( mean(combinedCope(:,baseline_time_indices,c,:),2) , [1 size(combinedCope,2) 1 1]);
                            %combinedCopeByState(:,:,c,:) = combinedCopeByState(:,:,c,:) - repmat( mean(combinedCopeByState(:,baseline_time_indices,c,:),2) , [1 size(combinedCopeByState,2) 1 1]);
                        end
                    end
                end
                
                % write back into the first level cope
                first_level_results.cope                = [combinedCope; magcope];
                first_level_results.stdcope             = [combinedStdCope; magstd];
                if 0,%NKglm>1,
                    first_level_results.cope_by_state       = [combinedCopeByState; magstatecope];
                    first_level_results.stdcope_by_state    = [combinedStdCopeByState; magstatestd];
                end;
                first_level_results.chanind             = first_level_results.chanind(1:size(first_level_results.cope,1));
                first_level_results.mask_indices_in_source_recon = first_level_results.mask_indices_in_source_recon(1:size(first_level_results.cope,1));
                
                % modify the sensorspace D object stored in first_level_results
                newcmblabels = newcmblabels';
                alllabels    = [newcmblabels; maglabels'];
                first_level_results.D_sensor_data = chanlabels(first_level_results.D_sensor_data,1:numel(alllabels),alllabels);
                nExcessChans = numel(first_level_results.D_sensor_data.chanlabels) - numel(alllabels);
                first_level_results.D_sensor_data = chanlabels(first_level_results.D_sensor_data,numel(alllabels)+1:numel(alllabels)+nExcessChans,[]);
                megcmblab = {'MEGCMB'};
                megcmblab = repmat(megcmblab,[numel(newcmblabels),1]);
                first_level_results.D_sensor_data = chantype(first_level_results.D_sensor_data,1:numel(newcmblabels),megcmblab);
                megmaglab = {'MEGMAG'};
                magmaglab = repmat(megmaglab,[numel(maglabels),1]);
                first_level_results.D_sensor_data = chantype(first_level_results.D_sensor_data,numel(newcmblabels)+1:numel(alllabels),magmaglab);
                otherlab = {'Other'};
                otherlab = repmat(otherlab,[nExcessChans,1]);
                first_level_results.D_sensor_data = chantype(first_level_results.D_sensor_data,numel(alllabels)+1:numel(alllabels)+nExcessChans,otherlab);
                
                first_level_results.D_sensor_data.save;
            case 'dont_combine'
                
        end % switch first_level.sensor_space_combine_planars
    end % if strcmp(source_recon_results.recon_method,'none')
    
    %% save results
    first_level_results.x=first_level.x;
    first_level_results.contrasts=first_level.contrast;
    first_level_results.source_recon_results_fname=source_recon_results_fname;
    first_level_results.first_level=first_level;
    first_level_results.source_recon=source_recon_results.source_recon;
    first_level_results.level=1;
    %first_level_results.source_recon_results_fname=source_recon_results.fname;
    first_level_results.name=first_level.name;
    first_level_results.recon_method=source_recon_results.recon_method;
    try
        first_level_results.session_name=source_recon_results.session_name;
    catch,
        first_level_results.session_name=source_recon_results.subject_name;
    end;
    first_level_results.artefact_map = artefact_map;
    
    % save results to disk
    if first_level.save_trialwise_data
        if ~isempty(first_level.trialwise_directory)
            trlwise_fname=fullfile(first_level.trialwise_directory, [first_level_results.session_name '_trlwise']);
        else
            first_level_results.fname=[ first_level_results.session_name '_' first_level_results.name ];
            trlwise_fname=fullfile(oat.source_recon.dirname, [first_level_results.session_name '_trlwise']);
        end
        disp(['Saving trialwise data in file ' trlwise_fname]);
        save(trlwise_fname, 'trlwise', '-v7.3');
        clear trlwise
    end % if first_level.save_data
    
    if first_level.doGLM
        first_level_results.fname=[ first_level_results.session_name '_' first_level_results.name ];
        disp(['Saving statistics in file ' oat.source_recon.dirname '/' first_level_results.fname]);
        osl_save_oat_results( oat, first_level_results );
        results_fnames{sub}=first_level_results.fname;
    end % if first_level.doGLM
    
    if (strcmp(source_recon_results.recon_method,'none') | isfield(first_level,'mni_coords')),
    else
        disp(['To create niftii files from this use a call to osl_save_nii_stats']);
    end;
end;


